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INTRODUCTION 


In  an  earlier  report  (Lyon  ^ » 1977),  a potentially  optimal  method  of 

recovering  deflections  of  the  vertical  from  RGSS  data  was  described.  This  report 
continues  that  work  — describing  the  implementation  of  the  method  and  estimates 
of  the  errors  associated  with  the  method.  In  the  first  section  of  the  report 
an  optimal  weighting  technique  is  derived.  This  technique  also  leads  directly 
to  a priori  error  estimates.  The  second  section  describes  the  results  from 
use  of  the  method  on  hypothetical  traverses.  From  these  data  it  appears  that 
the  optimal  method  can  indeed  lead  to  a significant  reduction  in  the  errors 
in  estimating  the  deflections  of  the  vertical.  A final  appendix  gives  instruc- 
tions for  the  use  of  the  associated  computer  program. 


J 


ERROR  SOURCES  AND  OPTIMAL  WEIGHTING 


There  are  three  sources  of  error  which  will  be  considered  in  this  report. 
These  are 

1)  correlated  gyroscope  errors 

2)  correlated  accelerometer  errors 

3)  errors  due  to  the  colocation  determination  of  the 
deflections  during  each  leg  - colocation  error. 

We  assume  that  the  gyro  and  accelerometer  errors  follow  a Langevin  equation 
(first  order  Markovian) 

^ + va  = A(t)  (1) 

where  v is  the  inverse  of  the  correlation  time  and  A(t)  is  Gaussian  white  noise. 
This  leads  to  a covariance  (Papoulis,  1965) 

< a(t^)  a (t^)  ^ exp(-v|t2  - t^|)  (2) 

2 

where  is  < o(t)  a (t)>,  the  variance  of  a.  Table  1 gives  the  values  for 

and  v used  for  this  report.  The  assumption  of  eq.  (1)  is  not  strictly 
true.  In  particular,  the  purpose  of  the  reduction  scheme  outlined  here  and 
in  the  previous  report  is  to  estimate  gyro  drift  rates.  The  deviations,  6a, 
from  this  estimate,  a,  say,  will  not  be  distributed  as  equation  (1).  However, 
it  is  clear  that  if  the  total  covariance  is  given  by  equation  (2),  then  t’r.a 
maximum  value  that  the  variance  about  the  mean  may  attain  is 

_ 2 

Var  (a-a)  < a (l-exp(-vT)) 

- 0 

where  T is  the  length  of  the  mission.  Similarly,  the  correlation  time  for 
values  about  the  mean,  t,  must  be  x < T.  Thus,  even  though  equation  (1) 
does  not  strictly  apply  to  variations  about  a,  we  will  assume  that  that 
variation  holds  and  the  6a(gyro)  are 

< 6a(tj^)  6a(t2)  > = a^  (1  - exp(-vT)  exp(-lt - t^  | /T)  (2a) 
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Table  1 


I 

[ 

Error 

! 

1 

[ 

Source 

Accelerometers 

Gyros 

I 


Parameters  Used  for  Deflection  Error  Estimates 


RMS  Value 


Correlation  Time 


10  microg's  (all  axes) 

2.5  X 10  ^ °/hr  (horizontal 

_’3  n 

2.0  X 10  ■ /hr  (vertical) 


AO  minutes 
2 hours 


I 


The  accelerometers  have  appreciable  white  noise  in  addition  to  the  correlated 
noise.  This  may  be  handled  approximately  by  increasing  the  accelerometer 
variances  and  decreasing  the  correlation  time.  The  values  in  Table  1 are 
based  on  the  data  given  by  Huddle  and  Maughmer  (1972)  suitably  modified  in 
accordance  with  the  discussion  given  above. 

The  colocation  errors  may  be  estimated  in  a straightforward  fashion.  We 
assume  that  the  actual  deflection  covariance  function  is  the  second  order 
Markovian  given  by  Kasper  (1971) . The  colocation  variance  is 

< (r  . - r.'^)  (r.  - r ®)  > = <rr.>  - < r ®r . ® > (3) 

iijm  ij  ij 

where  the  r notation  for  the  deflections  was  introduced  in  the  first  report, 
r.  is  either  a north  or  east  deflection  depending  on  whether  j is  odd  or  even. 
The  e superscript  denotes  the  estimated  value.  Equation  (3)  may  be  reduced  to 

<(r.-r.®)(r.-r.^)>  = <r.r.>-<r.r  ><r.r„><r.r  > ^ (4) 

1 i j 3 13  IK  3*’ 

where  the  tilde  denotes  deflections  belonging  to  the  basis  set  from  which  the 
others  are  estimated. 

The  basic  data  available  relate  to  u"  and  v*^,  the  north  and  east  velocity 
errors,  respectively,  at  the  end  of  the  n-th  leg  of  the  mission.  For  con- 
venience, introduce  the  notation 


n 


u 


n 


V 


(5) 


Then  the  basic  equation  for  the  error  velocities,  equation  (40)  of  the  original 
report  can  be  written  as 


n-1 


W.  - A.  p 
I £ic  K 


3=0 


S, 


(6) 
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where  A and  B are  matrices  defined  in  the  first  report,  p°  are  the  initial 
conditions  of  the  solution,  is  the  inhomogeneous  driving  term.  If  there 
are  no  errors  should  be  identically  zero.  If  errors  are  present  will, 

Xf  X- 

in  general,  be  non-zero.  To  estimate  the  errors,  we  square  equation  (6)  to 
obtain 


n-1 

, n.  2 , n n 0n2  , T,,,n  n o.  n j 

(F„)  = (W„  - A„„  p„)  + 2(W^  - A^^  M^)  ^ B 

J=o 


iK 


n-1  n-1 


+ Z E Bo".  B ” '"-i 

Horn  K o 


j=0  m=0 


(7) 


Ke  assume  no  errors  in  (w”  - and  that  the  expectation  value  of  the 

errors  in  ij;”,  - 0.  Taking  the  expectation  value  of  equation  (7)  gives 


2 n-1  n-1 

n 'S~'  \ ''  „ n „ n . ,m 

< F > / ^ . B„  . B.  < 6tii'^  > 

(i)  7^  f.Kj  Zorn  < o 

j=0  m=0 


(8) 


T i i 

where  = \ij''  (assumed)  - ii>-^  (true).  <F  > is,  of  course,  the  variance  of 

K K K » 

the  observed  data  point.  It  remains  only  to  evaluate  < 6iJ;^  >. 

Vc  use  the  ordering  of  i|;  of  the  first  report,  i.e.. 


= 


n 

gn 

0 

n 

a 


\-i 


I 


and,  then  6(1/  = 


-g<S.^"  + 6a" 


g6ri"  + 6a" 


6a 


6B 


6y 


(9) 
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respectively,  and  6a'^,  the  correlated  gyro  errors  for  Z,  N,  and  E 

axes,  respectively.  Neglecting  zero  cross-correlations,  the  error  covariance 
matrix  in  equation  (8)  becomes 


(Equation  10  on  following  page)  (10) 

Equations  (8)  and  (10)  then  give  an  estimate  of  the  errors  associated  with 

n n2 

the  measurement  of  W^,.  Furthermore  > is  the  optimal  weighting  for  the 

least  squares  solution  (Brownlee,  1962  ). 


P 
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ERROR  ESTIMATE  FOR  DERIVED  QUANTITIES 


We  need  now  to  derive  error  estimates  for  the  derived  quantities  in  the 
least  squares  solution.  We  rewrite  equation  (42)  of  the  original  report  as 


„n  „n  .no  „ , 

F.  =VJ„-A.  )J  ~D  h 
1 Z K 2.m  m 


(11) 


where  the  h are  the  quantities  in  which  we  are  interested,  i.e.,  the  calculated 
m 

deflections  and  gyro  drift  rates.  Change  the  notation  slightly  by  replacing 
the  double  index  (n,)!.)  by  j = 2(n-l)  + Z and  rewrite  equation  (11)  as 


F.  = a.  - D.  hm  (12) 

J J 

where  a.  = W.  - A.  The  matrix  equation  for  the  least  squares  solution  for 

J J J K K 

h is 


where 


and 


E..  h.  - b.  = 0 (13) 

J 1 


E 


(D  - D )(D  . - 5.)  W 

Ki  i KJ  1 K 


(14) 


D )W 

1 K 


(15) 


E 


D .W 

K1  Kj 


a 


N 


k:=1 


a W 

K K 


(16) 


2 

with  W^  = , the  optimal  weighting  discussed  in  the  previous  section. 

Writing  the  normal  equation  (13)  as  we  have  leads  to  a number  of  advantages 
(Brownlee,  1962  ).  The  solution  of  equation  (13)  is 

h.  = E.  “^b  (17) 

J J 1 i 
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and 


uu 


y 


(20) 


where  the  ' denotes  quantities  where  the  assumed  dependence  on  hy  has  been 
removed.  Thus,  if  we  wish  to  remove  the  gyro  drift  rates,  for  example,  from 
the  least  squares  solution  and  see  how  much  the  deflections  are  affectd,  it 
can  be  done  trivially. 

As  will  be  discussed  in  the  results  section,  we  follow  a somewhat 
different  procedure  in  eliminating  variables  from  the  least  squares  solution. 
Equation  (20)  h.olds  if  no  weighting  is  used,  or  if  the  weighting  is  unchanged 
after  removing  a variable  from  the  fit.  Since  we  remove  variables  that  do 
affect  the  weighting,  we  use  the  more  laborious  method  of  starting  from  scratch 
with  new  weights.  It  is  important  to  note,  however,  that  equation  (18)  still 
holds  and  provides  an  estimate  of  the  errors  involved  in  the  fit. 


.id 
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RESULTS  AND  DISCUSSION 


Two  hypothetical  traverses  were  used  to  find  estimated  errors  for  the 
outlined  reduction  method.  The  courses  are  sketched  in  Figure  1.  The  first 
traverse  is  a straight  line  to  the  northeast  covering  25  km.  The  second  is 
polygonal  - also  covering  25  km.  The  assumed  vehicle  speed  was  25  km/hr  - so 
that  total  travel  time  was  one  hour,  not  counting  stops.  Deflections  of  the 
vertical  were  determined  at  10  points  evenly  spaced  along  the  traverses.  The 
vehicle  was  assumed  to  stop  either  20  or  40  times  on  a mission.  This  made  the 
least  squares  system  well  overdetermined.  It  also  helped  produce  an  answer  to 
the  question  of  whether  fewer  or  more  stops  is  preferable.  Solutions  were 
generated  for  cases  including  all  the  gyro  drifts  as  fitted  variables,  includ- 
ing just  the  horizontal  axes,  and  including  none  of  the  gyros. 

Figure  2 shows  the  estimated  variance  in  the  north  velocity  channel  over 
the  course  of  a 20  stop  straight  line  traverse.  There  are  only  two  signifi- 
cant contributors  to  the  total  variance  - the  north  accelerometer  and  the  east 
axis  gyro.  Their  contributions  are  plotted  separately  in  Figure  2.  The  acce- 
lerometer error  is  more  or  less  constant  over  the  course  of  the  mission.  The 
gyro  drift  becomes  the  dominant  contribution  early  on  in  the  traverse  and  is 
constantly  increasing.  The  value  for  the  gyro  variance  used  for  Figure  2 is 
that  assuming  a constant  drift  rate  is  removed.  Without  the  removal  of  the 
average  drift,  the  gyro-related  variance  would  be  about  a factor  of  two  bigger. 
The  form  of  the  equivalent  curves  for  the  40  stop  traverse  is  nearly  identical. 
However,  the  individual  variances  are  about  half  what  they  are  for  the  20  stop 
case.  This  is  just  a reflection  of  the  fact  the  errors  at  individual  stops 
appear  to  accumulate  as  individual  random  events.  Half  the  time  then  implies 
half  the  accumulated  variance. 

The  results  for  the  polygonal  course  are,  once  again,  nearly  identical 
to  that  for  the  straight  line  traverse.  This  is  a direct  consequence  of  the 
fact  that  the  estimated  errors  introduced  from  collocation  are  negligible  in 
comparison  with  those  from  the  gyros  and  accelerometers.  According  to  the 
model  used  here  these  significant  error  sources  are  almost  independent  of 
the  direction  in  which  the  vehicle  travels.  Hence,  the  results  from  the 
polygonal  and  straight  traverses  are  almost  identical. 


-10- 


Figure  2.  Fst irmted  Variances  In  North  Channel  Over  25  kin  Traverse 


Figure  3 shows  the  estimated  errors  (standard  deviation)  for  the 
derived  north  deflection  of  the  vertical  using  20  stops  and  an  optimally 
weighted  solution.  The  first  point  to  note  is  the  terrible  performance  of 
the  fitted  solution  when  all  three  gyro  drifts  are  included.  Maximum  errors 
are  almost  200"  . This  behavior  occurs  only  in  the  north  direction.  In  the 
east  direction,  the  solution  is  as  well  behaved  as  the  other  two  curves  in 
Figure  3.  The  reason  for  this  strange  behavior  can  be  found  by  considering 
the  equations  for  the  east  gyro  error,  the  vertical  gyro  error,  and  the  error 
velocity.  Taking  those  equations  (7,  10,  and  12  from  the  first  report),  we 
find 

2 2 

r + r C-r  cos  d)  a + . . . (21) 

s £ s 

where  is  the  east  gyro  error,  r^  the  Schuler  frequency,  r the  terrestrial 
rotation  rate,  ^ the  north  deflection,  the  latitude,  and  a the  vertical  gyro 
drift  rate.  The  point  to  note  is  that  C and  a come  into  equation  (21)  in  the 
same  way.  Thus,  in  a least  squares  solution,  ^ and  a are  to  some  extent  inter- 
changeable. Since  there  is  no  vertical  channel  Information,  there  is  no  real 
way  to  separate  the  effects  of  a from  C-  The  east  deflection  is  well  behaved 
because  there  is  no  comparable  coupling  of  the  vertical  gyro  drift  rate  to  the 
east  deflection. 

A significant  improvement  is  made  by  removing  the  vertical  gyro  drift  from 
tile  solution,  as  can  be  seen  from  Figure  3.  The  results  are  still  somewhat 
puz; iing  as  the  estimated  errors  are  virtually  the  same  in  the  case  where  the 
horizontal  gyro  drift  rates  are  included  in  the  solution  as  when  they  are  not. 
This  is  in  spite  of  the  fact  that  the  assumed  variances  of  the  error  velocities 
is  about  a factor  two  smaller  when  gyro  drifts  are  included  in  the  solution. 
Unfortunately,  inspection  of  the  error  covariance  matrix,  eq.  (18-20),  shows 
that  ;;  and  y (the  east  gyro  drift  rate)  are  strongly  anti-correlated,  i.e. 
have  a large  negative  covariance  in  the  structure  of  the  least  square  solution. 
This  implies  that  the  situation  is  similar  to  that  discussed  with  respect  to 
the  vortical  gyro.  That  is,  with  the  given  information,  the  least  squares 
solution  has  difficulty  telling  the  difference  between  an  east  gyro  drift  rate 
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w 


and  a north  deflection  of  the  vertical.  This,  on  reflection,  should  not  he 
terribly  surprising.  We  can  write  an  equation  similar  to  equation  (21)  for 
the  error  velocity  in  the  north  direction.  Keeping  just  the  largest  terms, 
this  looks  like 


2 

r 


s 


u + 


Y 


(22) 


where  u is  the  north  error  velocity  and  y is  the  east  gyro  drift  rate.  The 
north  deflection,  C,  enters  equation  (22)  only  through  the  boundary  conditions. 
Another  way  of  looking  at  the  problem  is  that  equation  (22)  describes  a sinu- 
soidal variation  whose  phase  and  amplitude  depends  on  the  relative  sizes  of 
' and  Y-  In  effect,  the  phase  and  amplitude  must  both  bo  determined  by  one 
number  - the  error  velocity.  What  is  needed  to  help  is  information  about, 
for  example,  the  acceleration  error  at  a stop.  This  would  serve  to  disentangle 
the  two  quantities. 

Actually,  the  situation  is  not  nearly  so  bleak  as  has  been  painted.  The 
values  for  the  variation  of  the  gyro  drift  rate  about  the  mission  mean  arc 
quite  conservative.  The  actual  variance  could  easily  be  a factor  of  two  lower 
than  what  we  have  used.  In  this  case  the  solution  including  the  gyro  rates 
would  clearly  be  superior.  It  is  interesting  to  note  that  at  the  beginning 
of  the  mission  when  accelerometer  errors  dominate  the  variance,  the  two 
solutions  are  almost  identical.  This  implies  that  the  accelerometer  errors 
give  a limit  to  the  accuracy  of  the  recovery  of  the  deflection  in  the  neighbor- 
hood of  2-3"  for  the  20  stop  case  and  the  accelerometer  parameters  used. 

Somewhat  better  results  are  obtained  by  using  40  stops.  The  gain  is 
essentially  by  the  square  root  of  the  number  of  stops.  Thus,  a 40  stop  case 
gives  errors  about  a factor  V2  better  than  the  20  stop  case,  all  other  things 
being  equal . 

Since  the  major  sources  of  error  have  essentially  just  a time  dependence 
and  not  a position  or  velocity  dependence,  speeding  up  the  rate  of  traverse 
also  increases  tlic  accuracy.  The  increase  in  accuracy  is  roughly  proport ionn 1 
to  the  square  root  of  t!\o  velocity.  This,  of  course,  has  a limit  when  the 
vehicle  velocity  becomes  high  enough  to  make  the  neglect  of  velocity  dependent 
terms  in  the  error  propagation  equations  (first  report:  cq.  (1)  - (('))  serious. 
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So  far  the  results  discussed  have  dealt  with  weighted  least  squares 
solutions.  Figure  4 shows  a comparison  of  error  estimates  for  the  derived 
north  deflections.  Inherent  in  an  unweighted  solution  is  a single  assumed 
variance  for  all  the  error  velocities.  This  is  in  contrast  to  the  increasing 
- as  a function  of  time  - variances  in  the  weighted  solution.  It  is  not 
surprising,  then,  that  the  deflection  error  estimates  for  the  unweighted 
solution  tend  to  be  more  uniform  than  those  of  the  weighted  solution.  If 
the  error  model  used  is  reasonable  then  the  errors  derived  from  the  weighted 
solution  should  be  more  accurate.  In  practice,  the  derived  deflections  do 
not  seem  to  be  greatly  affected  by  the  choice  of  either  a weighted  or 
unweighted  solution.  Thus,  the  weighted  solution  appears  to  be  slightly 
preferable . 

The  results  presented  in  Figure  3 are  comparable  to  those  presented  by 
Huddle  (1973)  in  his  discussion  of  the  Position  and  Azimuth  Determining 
System  (PADS), 

We  have  argued  above  that  a factor  two  improvement  on  Figure  3 is 
easily  attainable  without  system  improvement.  This  would  be  superior  to 
the  PADS  results.  If  more  information  can  be  obtained  from  the  inertial 
system  - i.e.,  acceleration  errors  in  addition  to  velocity  errors  at  each 
of  the  stops  - the  accuracy  of  the  system  should  be  determined  by  the 
accuracy  of  the  accelerometers,  and  deflections  of  the  vertical  with 
accuracies  of  1 - 1.5"  should  be  attainable. 
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APPENDIX  1 


THE  OPTIMIZED  REDUCTION  PROGRAM 

A listing  of  the  FORTRAN  Program  used  to  determine  deflections,  and  error 
estimates  is  given  below.  The  input  data  are  described  in  the  comment  cards 
at  the  beginning  of  the  program  and  should  be  self-explanatory  with  two  excep- 
tions: 1)  all  input  data  are  cgs  and  angles  arc  in  radians,  and  2)  the  program 

is  set  up  to  handle  a traverse  with  known  deflections  at  the  start  and  stop. 

To  use  only  known  deflections  at  the  beginning  three  things  must  be  done.  First, 
add  a dummy  finishing  stop  to  the  data  V7ith  the  position  of  the  finish  equal 
to  the  start.  Second,  set  XIFIN  and  ETAFIN  equal  to  XI(|)  and  ETAij),  respectively. 
Third,  set  IDEFL  = 1. 

The  output  format  is  also  shown  below.  Solutions  are  given  for  cases  with 
all  gyro  rates,  horizontal  gyro  rates  only,  and  no  gyro  rates  in  turn.  Before 
the  actual  solution  the  estimated  variances  of  the  error  velocities  is  given 
both  in  total  and  from  the  individual  sources.  The  solved-for  quantities  are 
each  presented  with  error  estimates  (standard  deviations).  Finally,  for  each 
case,  the  actual  deviation  of  the  solution  from  the  data  is  given. 
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